#include "eigen.h"

Fwave::Fwave() {A.rl=0;A.im=0;B.rl=0;B.im=0;kv.rl=0;kv.im=0;mass=0;ptr=0;}
Fwave::~Fwave() {if(ptr!=NULL) delete ptr;}
Fwave::Fwave(Complex A_inp,Complex B_inp,Complex kv_inp,double mass_inp) {
  A = A_inp;
  B = B_inp;
  kv = kv_inp;
  mass = mass_inp;
}

Fwave::Fwave(Fwave & wave) {
  A   = wave.A;
  B   = wave.B;
  kv  = wave.kv;
  mass = wave.mass;
  ptr = 0;

}

Fwave Fwave::operator=(const Fwave & wave) {
  A   = wave.A;
  B   = wave.B;
  kv  = wave.kv;
  mass = wave.mass;
  ptr = 0;
  return *this;
}

void Fwave::is(Complex A_inp,Complex B_inp,Complex kv_inp,double mass_inp) {
  A = A_inp;
  B = B_inp;
  kv = kv_inp;
  mass = mass_inp;					
}

Complex Fwave::eval(double x) {
  Complex  coeff = cplx_sqrt((mass*MASS/PLANKUT)/kv);
  Complex ii(0,1);
  return A*coeff*cplx_exp(kv*x*ii)+B*coeff*cplx_exp(kv*ii*(-x));
}

void Fwave::show() {
  printf("A="); 
  A.show();
  printf(" B=");
  B.show();
  printf("\nkv=");
  kv.show();
  printf(" mass=%le\n",mass);
}


/* This routine finds number and position of the maximum of the function y(x).
	It also try to fit a lorenzian curve (sigma/( (x-energy)^2 + gamma^2 )
	for each found peaks.	
*/

Eigenvalue *  lorenzian( int num_y, double y[], double x[], int & num, int tail_inp, int Max_number_peaks_required)	{

	using namespace std;

  /* 	TAIL = Number of points from the peak on which we calculate 
	the coefficient gamma of the lorenzian curve.	*/

  int TAIL = 10, actual_number_peaks;
  int const MAX_PEAKS = 100;
  int k,j,i,count,lf,rh;
  double ratio,sum, param[MAX_PEAKS];
  Eigenvalue * z =0;
  string methodName ("Eigenvalue:lorenzian");

  if(Max_number_peaks_required>MAX_PEAKS) {
	  actual_number_peaks = MAX_PEAKS;
  }
  else {
	  actual_number_peaks = Max_number_peaks_required;
  }

  TAIL = (TAIL>=tail_inp)?TAIL:tail_inp;

  i = TAIL+1;
  count	= 0;

  while(i<num_y-TAIL-1)	{
    while( (y[i]<y[i-1] || y[i]<y[i+1]) && i<num_y-TAIL-1) i++;
    if(i<num_y-TAIL-1) {
      lf = i-1;
      rh = i+1;
      while( y[i]>y[lf] && (i-lf)<=TAIL ) lf--;
      if( i-lf>TAIL ) {
		while( y[i]>y[rh] && rh-i<=TAIL ) rh++;
		if(rh-i>TAIL && count<100) param[count++] = x[i];
      }
      i = rh;
    }
    i++; 
  }

  if(count > actual_number_peaks) {
	  count = actual_number_peaks;
  }
  num =	count;

  z = new Eigenvalue [num];
  if(z==NULL) {
	  cout << "OUT OF MEMORY LORENZIAN: Dynamic allocation failure for Eigenvalue " << methodName << endl;
  }

  i = TAIL + 1;
  for(k=0;k<num;k++) {
    z[k].energy = param[k];
    while(x[i]!=param[k]) i++;
    if(i>num_y-TAIL) {
      z[k].gamma = 0;
      z[k].sigma = 0;
    }
    else {
      j=0;
      sum	= 0.0;
      for(j=0;j<TAIL;j++) {
		ratio = y[i-j-1]/y[i-j];
		if(ratio!=1.0) {
			sum = sum+(ratio/(1.0-ratio))*(x[i-j-1]-x[i])*(x[i-j-1]-x[i])
			     -(1.0/(1.0-ratio))*(x[i-j]-x[i])*(x[i-j]-x[i]);
			ratio = y[i+j+1]/y[i+j];
			if(ratio!=1.0) { 
				sum = sum + (ratio/(1.0-ratio))*(x[i+j+1]-x[i])*(x[i+j+1]-x[i])
			     - (1.0/(1.0-ratio))*(x[i+j]-x[i])*(x[i+j]-x[i]);
			}
		}
		z[k].gamma = sqrt(sum/(2*TAIL));
		z[k].sigma = 0.5*y[i-3]*( (x[i-3]-x[i])*(x[i-3]-x[i])+z[k].gamma*z[k].gamma)
	           +0.5*y[i+3]*( (x[i+3]-x[i])*(x[i+3]-x[i])+z[k].gamma*z[k].gamma);
	  }
	}
  }
	
  return z;

}
